Stata: 面板数据模型-一文读懂

作者：游万海 || 连玉君 ||  (知乎 | 简书 | 码云)  

导言： 如下是连玉君老师上课的板书。你可以看出什么是 「固定效应」，什么是 「双向固定效应模型」，什么是 「POLS」 v.s. 「FE」 以及二者的差别。

所以，面板数据模型其实没有你想象的那么复杂！

常见的数据形式有时间序列数据( Time series data )，截面数据( Cross-sectional data )和面板数据( Panel data )。

从维度来看，时间序列数据和截面数据均为一维。面板数据可以看做为时间序列与截面混合数据，因此它是二维数据。数据形式如下：

世界是复杂的，所表现出来的行为特征也是复杂的，我们需要面板数据。

例如，欲研究影响企业利润的决定因素，我们认为企业规模 (截面维度)和技术进步(时间维度)是两个重要的因素。截面数据仅能研究企业规模对企业利润的影响程度，时间序列数据仅能研究技术进步对企业利润的影响，而面板数据同时考虑了截面和时间两个维度 (从哪个维度看都好看)，可以同时研究企业规模和技术进步对企业利润的影响。

正因为面板数据所具有的独特优势，许多模型从截面数据扩展到面板数据框架下。通过 findit panel data 命令可以发现目前Stata已有许多相关面板数据模型命令，包括(不限于)：

xtreg :普通面板数据模型，包括固定效应与随机效应        
xtabond/xtdpdsys/xtabond2/xtdpdqml/xtlsdvc：动态面板数据模型        
spxtregress/xsmle: 空间面板数据模型        
xthreg:面板门限模型        
xtqreg/qregpd/xtrifreg: 面板分位数模型        
xtunitroot: 面板单位根检验        
xtcointtest/ xtpedroni/xtwest: 面板协整检验        
sfpanel: 面板随机前沿模型        
xtpmg/xtmg：非平稳异质面板模型        

本文主要就普通静态面板数据模型进行介绍，包括模型形式设定、模型分类与选择及 Stata 程序实现等。

一. 模型形式设定

面板数据模型同时包含了截面和时间两个维度，设 i (i=1, ..., N) 表示截面 (个体)，t (t=1, ..., T) 表示时间，设定如下线性模型：

其中，

• 为 因变量，

• 为 自变量，

• 为模型误差项，是待估计参数，表示 的边际影响。

• 表示个体效应，表示那些不随时间改变的影响因素，如个人的消费习惯、企业文化和经营风格等；

• 表示时间效应，用于控制随时间改变因素的影响 (时间虚拟变量包括时间趋势项，时间趋势主要用于控制技术进步)，如广告的投放 (往往通过电视或广播，我们可以认为在特定的年份所有个体所接受的广告投放量相同)。

显然， 和 在多数情况下都是无法直接观测或难以量化的，因此也就无法进入模型。在截面分析中往往会引起遗漏变量的问题。

面板数据模型的主要用途之一就在于处理这些不可观测的个体效应或时间效应。当对所有的 i, 均相等时，模型退化为混合数据模型 ( Pooled OLS )，可直接用 reg y x 命令进行参数估计。

根据个体数N 和时期数 T的大小，通常可以将面板数据分为宏观面板和微观面板：宏观面板一般为 「大T小N」，微观面板一般为「小T大N」。依据 N、T大小不同，所采用的参数估计方法和分析中关注的重点也不尽相同。

二. 模型分类与选择

面板数据模型可以分为固定效应( Fixed effect model )和随机效应模型( Random effect model )。当 和相关，即，则该模型为固定效应模型；反之为随机效应模型。

两种模型的差异主要反映在对 “个体效应” 的处理上。

固定效应模型假设个体效应在组内是固定不变的，个体间的差异反映在每个个体都有一个特定的截距项上； 随机效应模型则假设所有的个体具有相同的截距项， 个体间的差异是随机的，这些差异主要反应在随机干扰项的设定上。

基于此，一种常见的观点认为， 当我们的样本来自一个较小的母体时，我们应该使用固定效应模型，而当样本来自一个很大的母体时， 应当采用随机效应模型。

然而，在具体的实例应用中，大母体和小母体并没有一个严格的界限，我们并不能明确地区分我们的样本来自一个较大母体还是较小的母体。因此，有些学者认为，区分固定效应模型和随机效应模型应当通过检验使用二者的假设条件是否满足。

下面我们讨论混合数据模型、固定效应模型和随机效应模型的选择。

2.1、固定效应的检验

固定效应的检验本质即检验个体间截距项的差异是否显著，即=0。根据假设检验原理，设定如下原假设

若结果拒绝原假设，则表明个体间截距项存在显著差异，模型中需要考虑固定效应。反之，混合 OLS 模型更为合适。通常可以利用 $F$ 统计量来检验上述假设是否成立：

其中：为固定效应模型的拟合优度系数(不受约束模型)，为混合数据模型的拟合优度系数(受约束模型)；N 和T 分别为截面与时期数；K 为解释变量个数。若原假设被拒绝，则说明个体效应显著，固定效应模型比混合数据模型更优。同理，可以构造相似的F统计量检验时期效应是否显著。

2.2、随机效应的检验

Breusch and Pagan (1980) 提出了基于面板随机效应模型残差的 LM统计量，构造如下原假设来检验随机效应：

相应的检验统计量LM为：

在原假设下，该统计量服从自由度为1的卡方分布。若拒绝原假设则表明存在随机效应。

2.3、固定效应还是随机效应？

通过检验说明个体效应 () 需要被纳入到模型中后，应该将 看成随机干扰项的一部分(随机效应模型)还是待估计参数 (固定效应模型)，下面介绍一些基本方法。

(1) Hausman 检验

从基本定义出发，可以通过通过检验个体效应与其它解释变量是否相关作为进行固定效应和随机效应模型筛选的依据。此时，我们可以采用 Hausman 检验。其基本思想是:在和其他解释变量不相关假定下，采用组内变换法估计固定效应模型和采用GLS法随机效应模型得到的参数估计都是无偏且一致的，只是前者不具有效性。若原假设不成立，则固定效应模型的参数估计仍然是一致的，但随机效应模型不一致。因此，在原假设下，二者的参数估计应该不会有显著的差异， 可以基于二者参数估计的差异构造统计检验量。

假设为固定效应模型的组合估计，为随机效应模型的 GLS 估计。在原假设成立下，有

根据方差公式

又因为，因此有

Hausman 检验基于如下 Wald 统计量

若拒绝原假设，表明个体效应 与解释变量相关，此时随机效应模型的结果不一致，应选择固定效应模型。

(2) 稳健Hausman检验Wooldridge (2002)

当不服从同方差假设时，传统的 Hausman 检验方法失效。Wooldridge (2002) 提出了一种稳健版的 Hausman 检验方法。建立如下辅助模型：

其中：为时变解释变量。当 RE 估计为完全有效估计时，利用 Wald 统计量做 检验所得结果应该渐近相等于标准的检验。当RE 估计为不是完全有效估计时，Wooldridge (2002) 提出在cluster-robust 标准误下做上述检验。

(3) 修正的 Hausman统计量

在固定效应模型与随机效应模型选择上，Hausman 统计量被广泛地应用于实证研究中。从上述看，该检验统计量渐近服从卡方分布，值应该为正数。然而，实际问题中计算出的统计值常出现负值的情况。针对出现负值这一现象，许多学者进行了研究，但并未形成一致的观点。

一种观点认为出现这样的情况主要是由小样本偏误引起，并建议此时应该解释为不能拒绝原假设，应选择随机效应模型 (如，Baltagi, 2008; Hsiao,  2003；Statacrop, 2009)。

另一种观点认为该统计量出现负值恰恰表明原假设不合理，此时应该选择固定效应模型。这些研究表明这种状况不仅仅出现在小样本情况下，在大样本情况下也时有发生 (Schreiber, 2008; Magazzini and Calzolarr, 2010)。如沈根祥 (2010) 在利用高频数据时也出现统计量为负值的情形。

连玉君等 (2014) 利用蒙特卡洛模拟方法得到内生性问题 (即解释变量与个体效应相关) 是导致统计量出现负值的主要原因。模拟分析表明，修正的 Hausman 统计量，以及过度识别检验方法能够很好地克服上述缺陷。

修正的 Hausman 统计量主要是对或进行调整。调整后的统计量为

或者为

其中：和分别为固定效应模型和随机效应模型下的均方根误差。

(4) 基于过度识别检验的 Wald 统计量

基于通常的 Hausman 统计量在存在异方差 (heteroskedastic) 情况下失效且当定义 cluster-robust 标准误时不再适用问题，Arellano (1993) 基于过度识别检验提出了 Wald 检验统计量解决这一问题。在条件同方差情况下，该检验统计量与通常的Hausman统计量渐近相等。此外，该统计量始终为正数。

如前所述，FE 估计和 RE 估计都需要满足一般意义上的外生性假设条件，即，而 RE 估计还要进一步满足面板特定的外生性假设条件，即。

我们可以将这个新增加的正交条件视为一个过度识别约束，以此来区分 RE 估计的前提假设是否合理。我们可以通过估计如下模型来构造 Wald 统计量

其中：，。和具有相似的定义。显然，在上式中，的 OLS 估计即为 RE 估计量，而的 OLS 估计即为之间的差异，即

利用 Wald 检验假设，所得统计量即为过度识别检验的 Wald 统计量。

(5) Mundlak’s (1978)  方法

在原假设成立情况下，估计量的有效性假设 (存在最小渐近方差) 是运用Hausman 检验的前提条件。然而，当误差项存在异方差或者序列相关时，这个条件往往不能够被满足。即使在这个条件满足情况下，该方法也可能存在小样本问题。 这里介绍另外一种方法，即 Mundlak’s(1978) 提出的一种检验方法。与通常的 Hausman检验不同，该方法在误差项不满足同方差和序列不相关情况下也是有效的。 设定如下线性模型：

Mundlak 方法的思想为检验和解释变量 是否存在相关。因此，建立如下关系式：

其中：是 的组内平均，是非时变的，且与自变量不相关的。

要保证 和解释变量 不相关，只需=0。根据以上式子，可以转化为检验如下方程的系数

因此，只需要回归这个方程，并检验是否成立。若拒绝原假设，则 和解释变量 存在相关，应选择固定效应模型。

(6) Bootstrap Hausman检验

传统的 Hausman 检验统计量可定义为

传统 Hausman 检验有效的前提条件是，在原假设为真的情况下，其中一个估计量为完全有效的。然而，实际应用中这个假设通常不被满足。特别地，当利用稳健标准误时，估计量通常非有效。

Bootstrap方法可以在估计量非有效的情况估计。假设重复进行 B 次抽样，可以得到 B 个和估计值，进而可得到 B 个估计值。可以利用下面式子进行估计

其中：。

三. Stata 实现

本部分以 Kleiber 和 Zeileis (2008) 的Grunfeld.dta数据集为例，说明运用面板数据模型的一般步骤。

3.1. 读取数据与面板数据设定

1. webuse grunfeld,clear //利用webuse从网络读取数据
2. list in 1/10          // 显示该数据集的前10行

1xtset company year,yearly //设置面板数据格式，利用 Stata 中`xt`开头的命令，必须用该命令进行设置。yearly表示年度数据，详细参考 `help xtset`

3.2. 模型检验与模型选择

本部分内容安排如下：

（1）个体效应和随机效应的联合显著性检验，以判别是否需要利用面板数据模型；

（2）若表明需用面板数据模型，利用Hausman统计量选择固定效应模型或随机效应模型更优；

（3）考虑到一般的Hausman检验在异方差和自相关情况下失效风险问题，对异方差，序列相关进行检验，以说明是否需要利用其它方法进行选择；

（4）针对一般的Hausman检验统计量可能为负值且对在异方差和序列相关情况不稳健问题，对稳健 Hausman 检验，修正的 Hausman统计量, 基于过度识别检验的Wald统计量法，Mundlak’s (1978)法，基于 bootstrap法的hausman检验等方法的Stata实现进行讲解。

（5）在选定固定效应模型或随机效应模型后，依据误差项结构（异方差，序列相关，截面相依）以及不同面板结构（「大T小N」，「大N小T」）， 介绍相应的参数估计命令。

（1）个体效应和随机效应的联合显著性检验

以invest为因变量，mvalue kstock为自变量，建立如下模型：

其中：和为待估系数。

利用Stata中 xtreg 可以方便实现面板固定效应模型与面板随机效应模型的估计。xtreg命令的语法如下：

1xtreg invest mvalue kstock,fe //fe表示固定效应;若同时包括时期虚拟变量，xtreg invest mvalue kstock i.year,fe，利用 testparm 检验                                   

1xtreg invest mvalue kstock,re //re表示随机效应

1xttest0  //检验随机效应是否显著，需要运行随机效应模型后使用

（2）Hausman检验

上述结果说明了有必要考虑个体效应和随机效应，接下来利用hausman 命令进行固定效应模型和随机效应模型的选择，主要步骤为：

• 步骤一：估计固定效应模型，存储估计结果；

• 步骤二：估计随机效应模型，存储估计结果;

• 步骤三：进行Hausman检验；

利用hausman 命令之前，有必要对其语法进行说明：

1. hausman name-consistent [name-efficient] [, options]

接下来进行hausman检验，

1xtreg invest mvalue kstock,fe
2est store fe_result
3xtreg invest mvalue kstock,re
4est store re_result
5hausman fe_result re_result

（3）异方差和序列相关检验

前文已经说明，当模型误差项存在序列相关或异方差时，此时经典的Hausman 检验不在适用，下面我们进行序列相关和异方差检验。

序列相关检验

先进行序列相关检验，在固定效应模型时可以利用命令xtserial，原假设为不存在序列相关。

1xtserial invest mvalue kstock

同样地，在随机效应时可以利用命令xttest1，原假设为不存在序列相关。

异方差检验

Greene (2000, p598) 提出一种修正的Wald统计量检验异方差，与标准的Wald统计量、LR和LM统计量不同，修正Wald检验同样适用于模型残差不服从 正态分布情况下。值得一提的是，在大N小T情况下，该方法的检验功效较低。该检验的原假设为同方差。

1xtreg invest mvalue kstock,fe
2xttest3

（4）模型选择其它方法

第一种：稳健 Hausman 检验**

目前 Stata 中没有相应的命令进行稳健 Hausman检验， 根据 2.3 中 (2) 部分公式，可以编写如下代码进行检验

 1webuse grunfeld, clear 
 2xtset company year
 3quiet xtreg invest mvalue kstock,re
 4scalar theta = e(theta)
 5global xlist2 invest mvalue kstock
 6sort company
 7foreach x of varlist $xlist2 {
 8     by company: egen mean`x' = mean(`x')
 9     generate md`x' = `x' - mean`x' 
10     generate red`x' = `x' - theta*mean`x'
11      }
12quiet reg redinvest redmvalue redkstock mdmvalue mdkstock, vce(cluster company)
13test mdmvalue mdkstock

第二种： 修正的 Hausman统计量

1xtreg invest mvalue kstock,fe
2est store fe_result
3xtreg invest mvalue kstock,re
4est store re_result
5hausman fe_result re_result,sigmamore

1hausman fe_result re_result,sigmaless

第三种：基于过度识别检验的Wald统计量

1 xtreg invest mvalue kstock, re cluster(company)
2 xtoverid 

运行后提示需要更高版本的ivreg2等命令，可以通过 net install ivreg2,from("http://fmwww.bc.edu/RePEc/bocode/i")进行更新。然后再运行

上述结果表明拒绝假设，应该选择固定效应模型。

第四种：Mundlak’s (1978)法

根据上文所述原理，可通过如下三个步骤实现该方法：

第一：计算解释变量均值

1local xlist "mvalue kstock" 
2foreach f of local xlist{
3bysort company: egen mean_`f' = mean(`f')
4}

第二步:估计包含均值的回归方程：

1xtreg invest mvalue kstock mean_mvalue mean_kstock,re vce(robust)
2est store Mundlak_result

第三步:利用test进行假设检验

1test mean_mvalue mean_kstock

结果如下

此外，也可以通过外部命令 mundlak 实现相同的系数估计，不过应该注意的是由于 mundlak不能得到稳健的标准误，得到的标准误和上述 手动运行方法不一致，所以test结果也就不一致。

1mundlak invest mvalue kstock,full
2test mean__mvalue mean__kstock 

第五种：基于 bootstrap法的hausman检验

由于存在序列相关和异方差，经典的hausman命令不再适用，下面使用基于bootstrap的hausman检验命令rhausman进行检验。

1xtreg invest mvalue kstock,fe
2est store fe_result
3xtreg invest mvalue kstock,re
4est store re_result
5rhausman fe_result re_result,reps(200) cluster

从检验结果可以发现，利用经典的hausman 和bootstraphausman均显示应该选择随机效应模型，而利用其他方法结果显示选择固定效应模型。

除了序列相关和异方差检验之外，截面相依检验也尤为重要。在固定效应模型中，可以利用命令xttest2进行检验，该方法是基于似不相关回归(SUR)进行 估计，所以一般要求截面数N比时期数T小；在随机效应模型中利用xtcsd进行检验，当然该命令也适用于固定效应模型。

（5）相关 Stata 命令推荐

依据误差项结构（异方差，序列相关，截面相依）以及不同面板结构（「大TT小NN」，「大NN小TT」）， 下文介绍相应的参数估计命令。

截面相依检验

1qui xtreg invest mvalue kstock, fe
2xttest2

1 qui xtreg invest mvalue kstock, re
2 xtcsd, pesaran

当误差项存在序列相关，异方差或截面相依时，依据形式不同，可以利用不同的方法和命令进行估计，详细可以参考 Hoechle (2007)。

几点说明

1. vce(robust) 和vce(cluster): 前者适用于异方差且观测值之间独立情况(heteroscedasticity-consistent standard errors)；后者 适用于异方差且允许观测值组内相关。例如cluster(group) 的含义是：假设干扰项在 group 之间不相关，而在 group 内部存在相关性。 若 group 代表行业类别，则表示行业间的公司所面临的随机干扰不相关，而行业内部不同公司间的干扰项存在相关性，或者是说，行业内的公司受到了一些共同的干扰因素。这部分内容将在后续的推文中详细介绍。

2. 固定效应模型与随机效应模型选择，学者们存在不同的观点。一些学者检验利用严格的统计检验选择，有些学者认为应该根据实际分析的需要进行选择，比如主要变量为不随时变的，那则必须采用随机效应模型。

3. 面板固定效应模型的估计除了可利用xtreg,fe进行估计外，也可以利用areg或者reg + dummy variables进行估计，注意这些方法的差异。

4. 上文中涉及到的一些命令，如xttest0, xttest1, xttest2, xttest3, xtserial, xtcsd, rhausman等需要下载安装。

4. 总结

虽然本文系统地介绍了静态面板数据模型的各种检验方法，但从现有的文献来看，实操层面的做法往往是单刀直入，甚至多少有些粗暴。

具体而言：

• 多数情况下 (90% 以上)，学者们都直接使用 FE，而 RE 则鲜有使用 (至少在公司金融和会计领域是如此)。

• 如果一定要在 FE 和 RE 之间进行筛选 (通常是为了应对审稿人)，建议采用假设较为宽松的 稳健 Hausman 检验 (help xtoverid) 或 bootstrap hausman 检验法 (help rhausman)。

• 在估计 FE 时，主流的做法是使用 「双向固定效应模型+聚类标准误」，即同时包含个体效应与时间效应的面板固定效应模型。对应的 Stata 命令为：xtreg y x1 x2 i.year, fe robust。注意：若仅关注系数估计值和其标准误，该命令等价于 xtreg y x1 x2 i.year, vce(cluster id) 以及 reg y x1 x2 i.id i.year, vce(cluster id)。换言之，xtreg, fe robust 中的 robust 选项本身就是在公司层面上聚类调整后的异方差稳健性标准误。

附录：文中所用 Stata dofiles

 1clear
 2webuse grunfeld,clear //利用webuse从网络读取数据
 3list in 1/10          // 显示该数据集的前10行
 4xtset company year,yearly //设置面板数据格式
 5xtreg invest mvalue kstock,fe //fe表示固定效应;若同时包括时期虚拟变量，xtreg invest mvalue kstock i.year,fe，利用 testparm 检验           
 6xtreg invest mvalue kstock,re //re表示随机效应
 7
 8xttest0  //检验随机效应是否显著，需要运行随机效应模型后使用
 9
10
11** 传统 hausman 检验
12xtreg invest mvalue kstock,fe
13est store fe_result
14xtreg invest mvalue kstock,re
15est store re_result
16hausman fe_result re_result
17
18xtserial invest mvalue kstock //序列相关检验,随机效应可以使用xttest1
19
20xtreg invest mvalue kstock,fe  
21xttest3                       //异方差检验
22
23
24** 稳健 hausman 检验方法
25
26quiet xtreg invest mvalue kstock,re
27scalar theta = e(theta)
28global xlist2 invest mvalue kstock
29sort company
30foreach x of varlist $xlist2 {
31     by company: egen mean`x' = mean(`x')
32     generate md`x' = `x' - mean`x' 
33     generate red`x' = `x' - theta*mean`x'
34      }
35quiet reg redinvest redmvalue redkstock mdmvalue mdkstock, vce(cluster company)
36test mdmvalue mdkstock
37
38
39**修正hausman检验方法
40xtreg invest mvalue kstock,fe
41est store fe_result
42xtreg invest mvalue kstock,re
43est store re_result
44hausman fe_result re_result,sigmamore
45hausman fe_result re_result,sigmaless
46
47
48
49** 基于过度识别检验法
50 xtreg invest mvalue kstock, re cluster(company)
51 xtoverid 
52
53
54** Mundlak’s (1978)法
55local xlist "mvalue kstock" 
56foreach f of local xlist{
57bysort company: egen mean_`f' = mean(`f')
58}
59xtreg invest mvalue kstock mean_mvalue mean_kstock,re vce(robust)
60est store Mundlak_result
61test mean_mvalue mean_kstock
62
63
64** 基于 bootstrap 法的 hausman 检验
65xtreg invest mvalue kstock,fe
66est store fe_result
67xtreg invest mvalue kstock,re
68est store re_result
69rhausman fe_result re_result,reps(200) cluster
70
71
72** 截面相依检验
73
74qui xtreg invest mvalue kstock, fe
75xttest2
76
77qui xtreg invest mvalue kstock, re
78xtcsd, pesaran

参考文献

1. 钟经樊和连玉君.计量分析与 STATA 应用,2010.

2. Hoechle D. Robust standard errors for panel regressions with cross–sectional dependence[J]. Stata Journal, 2007, 7(3):281- 312.

3. Breusch T S, Pagan A R. The Lagrange Multiplier Test and its Applications to Model Specification in Econometrics[J]. Review of Economic Studies, 1980, 47(1):239-253.

4. Mundlak, Y. On the pooling of time series and cross section data. Econometrica, 1978, 46:69-85.

5. Greene, W. 2000. Econometric Analysis.  Upper Saddle River, NJ: Prentice--Hall.

6. How can the standard errors with the vce(cluster clustvar) option be smaller than those without the vce(cluster clustvar) option? https://www.stata.com/support/faqs/statistics/standard-errors-and-vce-cluster-option/

7. https://blog.stata.com/2015/10/29/fixed-effects-or-random-effects-the-mundlak-approach/

8. Kleiber C, Zeileis A (2008). Applied Econometrics with R. Springer-Verlag, New York. ISBN978-0-387-77316-2, URL https://cran.r-project.org/package=AER.

9. Arellano, M. 1993. On the testing of correlated effects with panel data.  Journal of Econometrics, Vol. 59, Nos. 1-2, pp. 87-97. 10.Wooldridge, J.M. 2002. Econometric Analysis of Cross Section and Panel Data.  Cambridge, MA: MIT Press.

关于我们

• 【Stata 连享会(公众号：StataChina)】由中山大学连玉君老师团队创办，旨在定期与大家分享 Stata 应用的各种经验和技巧。

• 公众号推文同步发布于 CSDN-Stata连享会 、简书-Stata连享会 和 知乎-连玉君Stata专栏。可以在上述网站中搜索关键词Stata或Stata连享会后关注我们。

• 点击推文底部【阅读原文】可以查看推文中的链接并下载相关资料。

• Stata连享会 精彩推文1  || 精彩推文2

联系我们

• 欢迎赐稿： 欢迎将您的文章或笔记投稿至Stata连享会(公众号: StataChina)，我们会保留您的署名；录用稿件达五篇以上，即可免费获得 Stata 现场培训 (初级或高级选其一) 资格。

• 意见和资料： 欢迎您的宝贵意见，您也可以来信索取推文中提及的程序和数据。

• 招募英才： 欢迎加入我们的团队，一起学习 Stata。合作编辑或撰写稿件五篇以上，即可免费获得 Stata 现场培训 (初级或高级选其一) 资格。

• 联系邮件： [e-mail]

往期精彩推文

• Stata连享会推文列表1

• Stata连享会推文列表2

• Stata连享会 精彩推文1  || 精彩推文2

---